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The venerable 2D point-vortex model plays an important role as a simplified version of many 
disparate physical systems, including superfluids, Bose-Einstein condensates, certain plasma config- 
urations, and inviscid turbulence. This system is also a veritable mathematical playground, touching 
upon many different disciphnes from topology to dynamic systems theory. Point-vortex dynamics 
are described by a relatively simple system of nonlinear ODEs which can easily be integrated nu- 
merically using an appropriate adaptive time stepping method. As the separation between a pair of 
vortices relative to all other inter-vortex length scales decreases, however, the computational time 
required diverges. Accuracy is usually the most discouraging casualty when trying to account for 
such vortex motion, though the varying energy of this ostensibly Hamiltonian system is a potentially 
more serious problem. We solve these problems by a series of coordinate transformations: We first 
transform to action-angle coordinates, which, to lowest order, treat the close pair as a single vortex 
amongst all others with an internal degree of freedom. We next, and most importantly, apply Lie 
transform perturbation theory to remove the higher-order correction terms in succession. The over- 
all transformation drastically increases the numerical efficiency and ensures that the total energy 
remains constant to high accuracy. 
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I. INTRODUCTION 

A. Importance of Point- Vortex Dynamics 

Despite its relative simplicity, the 2D point-vortex model has appeared time and again in the description of a 
wide variety of physical systems. Part of this breadth of application is attributable to its long history starting with 
Hclmholtz in 1858 [T], who considered point-like distributions of vorticity imbedded in a 2D ideal, incompressible 
fluid. Indeed, the velocity field generated by the Hamiltonian motion of a collection of point vortices satisfies the 
Euler fluid equation [2]. This represents a remarkable computational and conceptual simplification: trading the 
infinite dimensional ideal fluid field equations for the finite dimensional coupled ordinary differential equations of 
the point-vortex model. It is completely natural then to ask: to what extent are various general phenomena in 2D 
hydrodynamics present in the point-vortex model? 

Onsager was one of the first to tackle this quiestion when he used ideas from equilibrium statistical mechanics 
to show the existence of negative temperature states [3]; states that correspond to large-scale, long-lived vortex 
structures, not unlike those that form in Earth's atmosphere. Further research in the statistical vein has resulted 
in a kinetic theory of point vortices ^4] , while consideration of structure formation dovetails nicely with the concept 
of Lagrangian coherent structures [5] [6j. There is also the view that point vortices are a useful toy model of 2D 
inviscid turbulence [7] [5] [5], an idea motivated by the non-integrable, i.e., chaotic, motion of four or more vortices. 
Closely related to this is the more rigorous concept of chaotic advection [Tn|, which can help explain the transport 
properties of a vortex dominated ideal fluid [llj . Clearly a wide variety of physical phenomena falls under the aegis 
of point-vortex dynamics. 

In addition to a diverse spectrum of phenomena in ideal fluids, the point-vortex model is also applicable to other, 
more exotic, physical systems. This includes 2D electron plasmas ^12j . Indeed, it was in the plasma physics com- 
munity where Lie perturbation theory was first used with great success [13] [M]. Bose-Einstein Condensates (BECs) 
also fall into this category. The Gross-Pitevskii equation, which governs the evolution of the BEC wave-function, 
can be re-expressed, via the Madelung transformation, as Euler's equations J.5 . Therefore, quantized vortex defects 
in a rotating BEC |16j |17j interact, to first approximation, as if they are point vortices. This also applies to other 
superfiuids such as He-II, where the vortex circulations are still quantized and the size of the vortex core is small 
enough to really warrant approximation by point vortices. Since supcrfluid turbulence is dictated by 3D quantized 
line vortices |18j , the point- vortex model can be considered a toy model of quantum turbulence as well. 

Aside from physical instantiations, the point-vortcx model should engender some intrinsic interest simply as an 
interesting mathematical entity. It has been described as a mathematical playground [1 , touching upon areas such 
as the theory of dynamic systems, ODEs, and Hamiltonian dynamics, whose appearance might be expected, as well 
as some ideas that at first seem to have no connection. For example, neither is it readily apparent that equilibrium 
configurations of point vortices can be connected to the roots of certain polynomials |19j . nor is it immediate that 
one can apply topology in the guise of Nielsen-Thurston and braid theory to describe fluid mixing [20 . A plethora 
of physical and mathematical considerations give weight to the notion that the point-vortex system is an important 
item of study despite its relative simplicity. 



B. Posing the Problem 

Orbits of the point-vortex dynamical system can be obtained easily by numerically integrating the set of coupled 
ODEs with an appropriate adaptive time-stepping method. When two like-signed vortices approach each other, 
they simply rotate about their center of circulation with an angular frequency that is inversely dependent on the 
separation length, squared. Therefore, in a system of many vortices, the closer a vortex pair is, the faster their 
angular movement will be compared to that of other vortices. To accommodate this motion, the integration method 
will reduce the time-steps between integrator function calls to maintain the prescribed tolerance. This prohibitively 
and unnecessarily slows down the integration of the system as a whole. One naive way of circumventing this problem 
is to adaptively change the tolerance, though this has the unpalatable consequence of decreasing the accuracy of 
the computed orbit. While one could possibly ignore this with appeals to the chaotic, non-integrable nature of the 
dynamics, there is another serious problem to consider. In most integrators the decreased accuracy will cause the 
vortices to systematically overshoot their ideal near-circular orbit. This has the effect of increasing the vortex-pair 
separation, and therefore decreasing the energy in this ostensibly Hamiltonian system. It turns out that one does not 
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have to choose either lowered accuracy and changing energy or long integration times. We develop a transformation 
in this paper which largely alleviates both of these problems, and does so in a physically intriguing way. 



Intuitively, a pair of like-signed vortices that are considered close when compared to all other inter- vortex distances 
- what we shall call a vortex "dimer" - will rotate about its center in a manner that appears unaffected by the 
presence of all other vortices. That is, the internal motion of the dimer, that of the constituent vortex pairs with 
respect to their center of circulation, is approximately the integrable motion of the pair on their own. Similarly, other 
vortices will interact with the dimer roughly as if it is a single vortex. This crude approximation has the useful feature 
that the now integrable dimer motion ceases to be the limiting factor in choosing the step size for numerical integration. 

Of course, we can not perfectly shoe-horn our problem into this picture. After some initial transformations to express 
the dimer in action-angle coordinates, we find correction terms which modify this picture. Fortunately these terms 
can be expanded in powers of a small factor, and therefore our system becomes amenable to perturbation methods. In 
particular, we use Lie transform perturbation theory, see [3T] [32] [23] [21], because the invariance of the symplectic 
structure of our phase space is manifest with these transforms. From here, we choose generating functions for the Lie 
transforms that will get rid of the correction terms successively at each order, while maintaining, as best we can, the 
integrable nature of the dimer's internal motion. When a correction term which we can not transform away does arise, 
we will see that it does not present any real computational problem. Indeed it will even have the interesting physical 
interpretation of an additional field that couples only to vortices which have an internal "spin" degree of freedom, 
i.e., the dimers. Overall we will have succeeded in converting the original N-body problem into an (N-l)-body problem. 

We will show that the overall transformation alleviates the small time-step issue, maintains a constant energy, 
and results in very accurate orbits. It should be noted that KAM theory and other superconvcrgent methods are 
not applicable here, seeing that they can not deal with the non-integrable system which arises after our first Lie 
transformation. Likewise it should be clear that our method is not connected with fast multipole methods, which 
increase the speed of integrating an N-body system by including only the most important pairwise interactions for 
each particle. We consider all N{N — l)/2 pairwise calculations, and therefore do not attempt to change how the 
algorithmic complexity scales with point-vortex number. We are more concerned with accuracy and how we can 
preserve it, while solving the specific small time-step issue. 



We consider a system of N point vortices in two spatial dimensions, each having position r j = {xj , yj) and circulation 
Tj . Their dynamics is described by the equations of motion. 



Here we have defined the * operator, so that if v — [v^, Vy) is a vector in two dimensions, then *v — {—Vy, v^), i.e., a 
rotation by |. These equations of motion can then be decomposed into coordinates, 



C. 



Outline of the Solution 



II. POINT- VORTEX DYNAMICS 



A. Equations of Motion 




(1) 




N 



(2) 




N 



(3) 



For reasons that will become clear, the 2iV-dimensional space whose coordinates are Vi,. .. ,ri^ will be called the 
phase space of the vortices. The state of all N vortices in the system is represented by a single point in phase space. 
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B. Hamiltonian Formulation 

The scalar function on phase space, 



N N 



i/(ri,...,r^.)^-^^^ln|r,-rfc|. (4) 



is called the Hamiltonian. It is easily verified that Eqs. ^ and ([s]) can be rewritten in terms of the Hamiltonian: 

dxj 1 dH 
~dt " Wj 
<Mj3_^_]_dH_ 

dt Tj dxj ^ ' 

Apart from the factors of l/Fj, these have the form of Hamilton's canonical equations of motion, where Xj and yj are 
canonically conjugate variables. 

Eqs. ([5) and ^ are rather unusual, in that most examples of Hamiltonian systems pair position coordinates with 
canonical y conjugate momentum coordinates. This system, by contrast, has pairs of canonically conjugate position 
coordinates. To pursue this identification, define the Poisson bracket of two functions 

r . ^ I ( OAdB dAdB\ 



J Tj \dxj dy, dyj dxj ^ 

It is readily verified that this bracket maps any two smooth (C°°) functions of ri, . . . , r^v to a third, is bilinear 

{aA + f3B,C}^a{A,C} + /3{B,C} (8) 
(where a and /3 are numbers), antisymmetric 

{B,A} = -{A,B}, (9) 

and obeys the Jacobi identity 

{{A, B},C} + {{C, A},B} + {{B, C},A} = 0. (10) 
In particular, we see that the phase-space coordinates themselves satisfy the bracket relations 

{xj,xk} ^ {yj,yk} ^0 (11) 
{xj,yk} ^ Sjk/^j- (12) 

Technically, because of the factor of l/Fj, the Poisson bracket defined by Eq. ([7| is noncanonical. It is possible to 
restore a canonical bracket by scaling the coordinates of the j''th vortex by factors of y^jT^f, swapping Xj and yj if Tj 
is negative f2i, but it seems easier to simply live with the noncanonical bracket that has emerged naturally from this 
analysis. All that is really required of a Poisson bracket is that it satisfy Eqs. (|8|, ^ and (10 1, and our noncanonical 
bracket does do so. 

In terms of the Poisson bracket, the canonical equations of motion may be written 



dxj 
~dt 

dt 



^{x„H} (13) 



Indeed, the rate of change of any function on phase space, A{ri, . . . , rjv), is given by 
dA f OA dxj OA dy, \ ^ 1 f dA dH dA dH 



3 ■' 3 



dt ^\dxj dt dyj dt ) ^T.j \dxjdyj dy^ dxj ) ^ ' ' ^ > 
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Conservation of energy is then established by noting that 

^ = {i/,i/} = 0, (16) 

which follows from the antisymmetry of the bracket. That is, the Hamiltonian is constant, and numerically equal to 
the vortices' energy, 

H{r,,...,rN)=E, (17) 
constraining the dynamics to lie on a surface of codimension one in phase space. 

III. A VORTEX DIMER 

A. Hamiltonian for a Vortex Dimer 

Returning to Eq. Q for the Hamiltonian of N vortices in an infinite domain, we suppose that two like-signed 
vortices are very close to one another. We further suppose that these two are the ones labeled m and n, and we refer 
to this pair as a "vortex dimer." We rewrite the Hamiltonian as follows: 

r r ^ ^ Y Y 

i?(ri,...,r„i,...,r„,...,rjv) = In |r,„ - r„| - ^ ^ --^\n\rj - rk\ 

E ^lnK.-r,|- ^^ln|r„-r,|. (18) 

The first term on the right is the Hamiltonian of interaction between vortices m and n within the dimer. The second 
term is the Hamiltonian of interaction amongst all of the other vortices in the system. The third and fourth terms 
describe the interaction of the dimer with all the other vortices in the system; in particular, the third term describes 
the interaction between all of the other vortices and vortex m, and the fourth term describes the interaction between 
all of the other vortices and vortex n. 

B. Transformation to Center of Circulation and Relative Coordinates 

We introduce the center of circulation of the vortex dimer 

^= r +r — ' 

and its relative displacement 

r = rn- Tm- (20) 

We transform coordinates to eliminate and r„ in favor of i? = (X, Y) and r = (x, y). The inverse transformation 
is 



r„r = R- ^ ^"^ r (21) 
rn = R+ ^ r. (22) 



The Poisson bracket relations among the new coordinates may readily calculated to be: 



= (23) 



{x,y}^^, (24) 
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and 



{x,y} = {2/,x} = o, 



where the total circulation of the vortex dimer is 



and the reduced circulation of the vortex dimer is 



r + r 

It follows that the Poisson bracket in the new coordinates is 

1 fdAdB dAdB\ 1 fdAdB dAdB 



^ 1 / dA dB dA dB 



(25) 
(26) 

(27) 
(28) 



We note that X and Y comprise a canonically conjugate pair, as do x and y. 



It remains to write the Hamiltonian in the new coordinates. Using Eqs. (21 1 and (22), this is straightforward, and 
we find 



N N 



r.r 



k 



H = — In r — y , 

27r ' ' ^ ^ 47r 



- y i^ln 
27r 



R Tj —r 



ln|rj - rfcl 



27r 



In 



^ In I r I — — - In I r , 



^ r r 

- y i^ln 
^ 47r 

^ r V 

- J2 



i m. 



r_Rrr 

27r 

AT 



47r 



In Irl 



2^(l?-r,) 



TV 



E 



TV 



r2 



rr 



In I r , 



^ ^™^ln| 1-2,^'^ ^^'"'^ 



rk\ 



j^m,n 



27r 



N 



4tt 



^ ^"^Mn(l + 2r ^""'^^ 



47r 



\R-r, 



(29) 



As we will be referring to the terms of this Hamiltonian often, it makes sense to label them. The first three terms 
are collectively denoted Hq. The first term on its own is Hqi, while Hq2 corresponds to the second and third terms 
together. The Hamiltonian of interaction between the two vortices within the dimer is i?oi- The first term of Hq2 is the 
interaction Hamiltonian of all non-dimer vortices in the system, and the second term of Hq2 is the principal interaction 
Hamiltonian between the dimer and all of the other vortices. Altogether Hq makes the crude approximation that 
the dimer, while having internal degrees of freedom, is merely a vortex of circulation Tn at position R, interacting 
normally with the N-2 other vortices. The rest of the Hamiltonian constitutes correction terms to this picture. 



IV. ORDERING 



A. Ordering the Hamiltonian 



We would like to develop an ordering scheme whereby the fastest and most important motion is the oscillation of 
the relative displacement vector of the vortex dimer, described by the first term (Hqi) in Eq. (29). We introduce a 
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formal ordering parameter, e, where the small quantities, |r|/|i2 — rj\ and iT-j/jri 



are of order e or smaller. 



Eq. ( 29 ) is thus rewritten 



N N N 



27r 

N 

E 

N 



r r 

47r 



In 1 - 2e 



I R - 



r r ■ 

47r 



(R-Vj) ■ r 



(30) 



One should immediately note that in this scheme, all of Hq is of order zero in e, and therefore the internal dimer 
motion alone does not yet enjoy the distinction of being the unperturbed motion. We will soon remedy this, when 
we also order the Poisson bracket, or equivalently the Poisson tensor, in e. The Hamiltonian ordering does, however, 
indicate that the correction terms are of order and higher, though it is not yet manifest. 



B. Poisson Tensor 



Lie transform perturbation theory will require that we use the Lie derivative, see Eq.(64|. To ensure that the Lie 



derivative will deal with the e ordering correctly, we will need to order the Poisson bracket as well as the Hamiltonian. 
For the sake of transparency we adopt the Poisson tensor framework of symplectic geometry instead of the Poisson 
bracket formalism. These are both interchangeable, but ordering, and indeed perturbing, the Poisson tensor is 
conceptually much simpler. We can define can define the Poisson tensor P in terms of the Poisson bracket: 



(31) 

This Poisson tensor has constant components, see Eq.(|35| below. With it we can rewrite the equations of motion as 



which we will usually refer to in the following equivalent shorthand notation 

z = VdH. 



(32) 



(33) 



C. Ordering the Poisson Tensor 

We have two goals that we wish to accomplish with the ordering of the Poisson tensor. The first goal, alluded to in 
the previous section, deals with how the Lie derivative handles the ordering. This primary concern will be explained 
later, after we have introduced the Lie derivative machinery. However, we will state here the ordering of the Poisson 
tensor that arises from this consideration: 



where the (Q)-tensors Pq and P2, expressed as 2N x 2N matrices, are: 



A 

-A 



B 

-B 



where A ^ 
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and 



lo 



(34) 



0\ 



(35) 



Our second goal is to recognize that the fast internal movement of the dimer is the most important motion, and 
therefore elevate it to the status of unperturbed motion in our eventual Lie transform perturbation theory. The 
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Poisson tensor ordering of Eq.(34) achieves this. That is, the zero-order term of i = VdH, when expanded in e, is the 
motion of the dimer alone: 

i = (Po + e^Ps) d {Hn + 0{e^)) = Porfi^oi + e^P2dHo2 + O(e'), or (36) 

While this justifies our consideration of the internal dimer motion as the unperturbed motion, one should take 
care in interpreting the rest of the above dynamic expansion. In particular, the ordering is absolutely correct only 
when comparing the relative size of different contributions to the evolution of an individual coordinate, and not when 
comparing a single contribution to the evolution of one coordinate with that of another coordinate. This is a minor 
issue, and certainly will not affect the development of our overall Lie transform method. 

V. THE UNPERTURBED PROBLEM 
A. Action-Angle Variables 



The Hamiltonian of the unperturbed problem, (i?oi) of Eq. (30) is 



Ho^ = -^ln\r\. (38) 

Note that this depends only on the relative displacement of the dimer r, so x and y are the only coordinates that 
vary along the unperturbed orbits. Their equations of motion are 



1 OHq ^ Tn V 

Tr dy 2tt + y"^ 

I dHp ^ ^ Tr X 

Tr dx 2tt x^ + y^ 



(39) 
(40) 



It immediately follows that 

d 



^ {x^ + y2) = 2xx + 2yy - 0, (41) 



so that |rp = + is a constant of the unperturbed motion. 

We define the action variable J and the angle variable 6 of the unperturbed motion, 

J=^\r\' = ^{x' + y^), e = '^rg{y + tx). (42) 

We transform coordinates again, this time to eliminate x and y in favor of 9 and J. The Poisson bracket relations 
among the new coordinates may be calculated as follows, 

1 fdOdJ dedJ\ If y -x \ J_ 

It follows that 6 and J comprise a conjugate pair of coordinates, and that the Poisson tensor and its ordering remain 
unchanged if the coordinates are written, 

z = {e,X,x„--- ,J,Y,y„---). (44) 
We next aim to write the Hamiltonian in the new coordinates. Toward this end, we define the angles 

0,^arg[{Y-y,) + i{X-x,)], (45) 

so that 

[R - rj) • r = V2J \R - rj\ cos {9 - 6j) . (46) 
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We also make use of the straightforwardly derived Fourier expansion 



In (1 =F 22; cos 6* + 



cos 



z . 



Then, after some algebra, the last two terms of Eq. ( 30 1 transform to give us 

N 



1, N N 



47r 



+ (-i)^r 



e-1 



N 

E 



r r 



An 



£=2 

We have thus succeeded in writing the Hamiltonian as 

C30 

where Hq = -ffoi + -^02, 



2J \ 



2n 



\R- 



cos [e {9 - Oj 



Hoi 



An 

N 



In (2 J) 

N 

E E 



An 



N 



ln|rj-rfc|- ^ 



and for i > 2, 



rl-^ + (-i)^r 



^ r r 

An 



2J 



\R 



cos 



(47) 



(48) 

(49) 

(50) 
(51) 

(52) 



The rather simple form of the higher-order terms Hi is encouraging. To the extent that the unperturbed motion of 
r is oscillatory, these terms will be straightforward to average and integrate over unperturbed orbits. 

It is very important to note that and J are action-angle variables for the unperturbed problem with Hamiltonian 
Hqi. They are not action-angle variables for the full Hamiltonian H. Our strategy now will be to use perturbation 
theory to find a near-identity canonical transformation so that 9 and J are action-angle variables to higher order in e. 



B. The Unperturbed Motion 

The Hamiltonian Hqi is independent of all coordinates other than J, and therefore the coordinates R, J, and rj 
for j 7^ m, n are constant along unperturbed orbits. Only 9 varies along unperturbed orbits according to 

9 = FodHoi=V'o'^--^- (53) 

It follows that 

9{t) = 6*0 - nt, (54) 
where is a constant of integration, and where we have defined the rotation frequency 

O.i^. (55) 

This means that averages of a phase function A along unperturbed orbits are simply averages over the angle variable 
9. These are accomplished with the operator 



{A) = ^ I A{9) d9. (56) 



r2TT 
2^ Jo 

The oscillatory part of a phase function A is then denoted 

A = A^{A). (57) 
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^^01 


= 




= 


He 


^Hi 



C. Averaging the Hamiltonian Over Unperturbed Orbits 

We are soon going to need the average of our Hamiltonian over unperturbed orbits, so we compute it here. Because 
Hqi and Hq2 are independent of 9, they are unchanged by averaging over 0; that is 

(ffoi) - Hoi, (58) 
{H02) - H02. (59) 
For i>2, since (cos [£ {9 - 0j)]) = 0, it follows that 

(Hi) = 0. (60) 

The oscillatory parts of the Hamiltonian are then 

(61) 
(62) 
(63) 

VI. LIE TRANSFORM PERTURBATION THEORY 
A. Overview 

At this point, our phase-space coordinates are z = {9, X,Xj, ■ ■ ■ , J,Y,yj, ■ ■ where j ^ m, n. We introduce the 
Lie derivative £i}, with respect to the vector field v. It acts on scalar fields / (z) and on other vector fields w (z) in 
the following manner: 

£vf^v'f,:^v'^^, (64) 

i£vw f = v'w^, - w'v% (65) 

We can take the Lie derivative of any tensor by contracting this tensor with the requisite number of arbitrary 
1-forms and vector fields, applying the Lie derivative to the resultant scalar, and then applying the Leibnitz rule. We 
will need the Lie derivative of a Poisson tensor 

i£vVy'' = V'J'v' - P*'=i;f, - r'v% (66) 



but the first term in Eq. (66) is equal to zero, since our Poisson tensor has constant components. 



The method of Lie transformations usually proceeds by picking a generating function, g, which gives a vector field 

= V"'9,k (67) 
and noting that the following coordinate transformation is always canonical 

z = e-xp {+£i}) z. (68) 
That is, the Lie derivative of P disappears under such a generating function, thereby preserving the symplectic 



structure. While this choice, Eq. (67), overdetermines the generating vector field, we will use it until a little more 



subtlety is need to address fourth-order corrections. This transformation results in a new Hamiltonian 

H ^eyiY>{-£v)H. (69) 

If V is of order e-' , the transformation will affect the Hamiltonian only at order and higher. 

Our strategy will be to determine successive vector generators v,;, starting at second order, that preserve the Poisson 
tensor and eliminate the Hamiltonian at all higher orders. To leave the largest variety of near identity transformations 
at our disposal, we choose an infinite product form, in the manner of Dragt and Finn |21j : 

z = exp (+£^£2) exp {+e'£z) exp {+€^£4) ■■■z. (70) 

Here, = £vn- This results in a new Hamiltonian and Poisson tensor: 

77 = • • • exp {-e^£A) exp {~e'£z) exp {-€^£2) H, (71) 

P = • • • exp {-e^£A) exp {~e'£z) exp {-€^£2) P- (72) 
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B. Ordering 



Poisson Tensor Ordering Redux 



Now that we have the Lie derivative machinery at our disposal, we can justify our use Eq.p4[) for the initial Poisson 
ordering. Consider the Lie derivative of a scalar /, Eq.(64|, with respect to a vector field, Eq.(67), generated by the 
scalar g: 



£vf 



ik dg df 



(73) 



Because of the symplectic structure of the Poisson tensor, we can partition the terms on the right hand side into two 
groups: those that have derivatives with respect to variables (0, J) , and those that have derivatives with respect to 
variables (X, , • • • , 1", ?/j , • • • ) . For the Hamiltonians and scalar generators we are considering, the magnitude of the 
{9, J) terms in the Lie derivative will be reduced by roughly a factor of J when compared to the magnitude of g and / 
together. Similarly, the magnitude of the (X, Xj, - ■ ■ ,Y,yj,- ■ ■) terms in the Lie derivative will be reduced by roughly 
a factor of |-R — rj\ or jr, — rj\ . Therefore the ratio of the {X, Xj, - ■ ■ ,Y,yj,- ■ ■) terms in the Lie derivative to the 
{9, J) terms will be of order e^. This relative ordering can be formally achieved by ordering the initial Poisson tensor 
in the manner given by Eq. ( 34 1 : 



£vf 



dz^ dzi 



ijfc dg df 



dz^ dz^ 



(74) 



2. Lie Transform Ordering 



Expanding each of the exponentials of Eq.(71 & 72 1 to fourth order and putting in the explicit initial ordering of 
the Hamiltonian and Poisson tensor gives: 



i? = (1 - e^A) (1 - e^£,) [l - ^£2 + e'\£i) (h, + j 
P = (1 - e^£,) (1 - e^£,) [l - e^£, + e^l^i) (Pp + e^V,) 



(75) 



(76) 



We suppose that both the final Hamiltonian and Poisson tensor are also ordered in e. Then, expansion of the above 
in powers of e and matching terms yields the sequence of equations 



Hf) — Ho 

H2 = H2 — £2Ho 
H3 = -ff 3 — £3Ho 



(77) 
(78) 
(79) 
(80) 

(81) 



and 



= 
= ^£^ 



£2 



-£4^0 — £2^2 + 2'^2^o 



(82) 
(83) 
(84) 
(85) 

(86) 



After applying the Lie transforms, we would like the equations of motion to consist of the motion of the dimer and 
the motion of all the other vortices with the dimer. That is, z = VodHo + P2(ii?o- We can ensure this by requiring 
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that order by order, starting at second order, the transformed Poisson tensor remains unchanged and the transformed 
Hamiltonian disappears. Furthermore, secular perturbation theory forbids us from including any averaged, Eq. (56 1, 



contributions to our transformed Hamiltonians in the calculation of the vector generators. This means that we will 
need to check that {H^) = at each order. These considerations lead to the following conditions: 



H2 — =f 2^0 — 
-^3 — ^^Ha = 



and 



£2^0 = 

i^gPo = 



£7 



1 



2J^0 



0. 



(87) 
(88) 
(89) 

(90) 
(91) 

(92) 



3. Orders Two and Three 



To maintain the second-order motion that results from the interaction of the non-dimer vortices, as well as their 
interaction with the dimer as though it were a single particle of circulation F^j, we must demand that the conditions 
in Eq. (87) & (90 1 are met. The second of these conditions is trivially fulfilled if 



(93) 



Thus the vector generator at second order, which is now completely determined by a scalar generator 52, has only 
two non-zero components. 



V2 



I 1 dg2 
'Tr dJ 



1 dg2 
Tr 89 



0,- 



This allows us to write the condition, Eq. (87), on our Hamiltonian as: 



H2 



Tfl dg2 

An J ae 



0. 



(94) 



(95) 



If we demand that g2 be single- valued to preclude secular terms and consider Eq. ( 60 ) , it is immediately apparent 
that (H2) = {H2) - w) = 0- ™ust solve Eq. ^ for the generator 
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AttJ 



H2 [9') d6'. 



From Eq. ( 52 ) , this yields the generator 
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J' 



N 



sin {29 - 29-1 



(96) 



(97) 



At this point it would be ideal to be able to renormalize the unperturbed orbits by treating all of Hq as if it were 
the new unperturbed problem, in the manner of a superconvergent Lie transformation. However, this problem is 
already non-integrable and therefore not amenable to this treatment. So, we proceed to third order with the same 
unperturbed motion. Identical considerations lead to the following third-order generator 
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V2Tr (F„ — F,„) 



F2 
^ R 



J5 



N 

E 



F. 



sin(36i- 36ij ) 
\R-rA' 



(98) 



The generator 172 provides the most important correction term present. With it, the transformation given by 
Eq. ( 70 ) takes us to a set of coordinates in which the dimer may be approximated as a single vortex of circulation 



Tfj and position R. Whereas this approximation was valid only to order e in our original coordinates, in the new 
coordinates it is accurate to order - that is, the first corrections to it are order e'^. When the generator 53 is 
included the validity of this approximation is extended to order e'^. We next refine this transformation to push the 
corrections to still higher order. 
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4- Order Four 



Using £2^0 = and H2 — =£2^0 — 0, the conditions at fourth order, Eq. (89 1 & (92), can be written 

1- 



H4 — £iHi^ — -£2H2 — 



(99) 



and 



£4P0 + £2P2 = 0. 



This latter equation is satisfied with the choice 



Our vector generator at fourth order now has the form 



,J^dg4 J^dg2 1 dg2 



1 dgi J^dg2 1 9.92 
' de ' Tr dX ' Tj dxj 



(100) 



(101) 



(102) 



All but two of these terms are already determined from g2- We will fix 54 by considering Eq. (99). Once again 
we demand that 54 is single- valued and check for any averaged contributions to -ff4. Unlike at second and third 
order, we encounter an averaged contribution to H4, due to ■£2-f^2, that secular perturbation theory forbids us from 
transforming away. Before we consider how this averaged term affects our perturbation theory, we can derive 34 from 
the purely oscillatory parts of each term: 



54 



Ha - Pf 52,fci?0j- - ^Pf'.92.fci/2 J- 



d0'. 



From this the fourth-order generator is straightforwardly, though laboriously, derived: 



(103) 



N 



9i 



4r| 



p3 



^ R 



sin (461 - AO J 



8(r; 



sin(26l-26lj) 



N 



sin {20 -39 j + 0k) 



\R-r,\^\R-rk\ 
where 0jk = arg [{yj - yu) -f i {xj - Xk)]. 



sin [20 - 30 j + 0jk) sin (461 - 261^ - 26ife) ' 
\R-rA'\r,-rk\ " \R - rA' \R - r^f. 



(104) 



5. Averaged Contribution to Hi 

At fourth order, we must consider the effect of the non-oscillatory term due to £2^^2- This gives us an averaged 
contribution to the Hamilton, 

, _ 3T,..P ^ ^ cos (20, -20,) 

""'^^ j^„i,nk^m.,n I""" \"- 

which changes the dynamics 



(106) 



Only the angle variable of the dimer, 0, is affected at fourth order. In particular, J remains a constant of the motion 
and 

= -n + e^^ E E r.r^ ^os(2^,- - 2g,) 

27rr« .f- f- ' \R-rA^\R-rk\^ 
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Thus the rotational frequency of the dimer decreases at fourth order. Though the extent to which it does decrease 



depends on the configuration of other vortices. To help interpret Eq. (107), we can write the unit vector from vortex 



k to the dimer as Rh 



R-r, 
\R-r, 



and use the * operator from Eq. ( 1 ) , equivalent to rotating the vector by f , which 



gives the following for the important part of the frequency correction: 



cos {20 J - 29 k) 
\R-r,\'\R-ru\' 



Ry[R, 



Rl 



Rk 



\R^ 



IR'Vkl 



R 



(108) 



Though it is not obvious, any configuration of any number of extra vortices will always decrease or maintain the 
dimer's rotational frequency. For a single extra vortex, with either positive or negative circulation, this is apparent. 
The complications arise when considering the three-body contributions to Eq.(lG7). We can see, from Eq.(108), that 



two extra, like-signed vortices which are co-linear with the dimer will further decrease the rotation rate. If, however, 
they are both an equal distance from the dimer and are at right angles to each other with respect to the dimer, there 
will be zero net change in the dimer's rotation rate. The results of these two scenarios are interchanged when the 
two extra vortices have opposite signed circulations. When there are more than two extra vortices, it becomes more 
difficult to argue that the change in the dimer's rotation rate is negative definite, though we have not seen a counter 
example in numerical simulations. 



As a physics analogy, we consider 9 to be an internal variable of the dimer, much like spin. Up until now the 
rest of the system has been transparent to the spin of the dimer. With the addition of the fourth-order averaged 
Hamiltonian, we can view the spin as coupling to a new field; a field generated by the vortices, but coupled to only 
by vortices with non-zero spin, i.e., the vortex dimers. This new field gives the dimer a position-dependent effective 
spin, which is smaller than the bare spin {D,). This interesting view remains valid through at least fifth order in our 
calculations. 



C. Transformed Variables 



To fourth order, the transformation given in Eq. 70 is written as 



(109) 



2 "p2 ~ Ej 1 and (3j = (Tn -f E^). To get the 
forward transformation, simply use the upper operator of ± or ^; to get the backward transformation, use the lower 
of ± or =F and replace all transformed coordinates z with the original coordinates z and vice-versa. 



N 



sin {29 - 29 j) 



^ 3 10^2 (E„ - E™) 3 ^ ^sin{39-W,) 
^ ^ ^ 2^ ^J^TT, — 1713^ 



E^ 

^ R 



\R-r, 



N 



3 \ sin (461- 461 
±E, + -a, > ^ ^' 



R - r. 



sin (26*- 261, ) 



N I 



sin {29 - 39 J + 9k) ^sin {29 - 39 j + 9jk) 
\R-r,f\R-rk\ ~ |H - r,f |r, - r,.| 

3^^ sin(46'-26'j -26*^)^ 



4/ \R-rA^\R-rk\' 



(110) 
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N 



cos (26* - 20 J 



N 



3 4V2(r„-r„0^| j2 p.cos(30-30j) 



^ R 



N 



T ^^^^^ E 



cos (461 - 46ij 



\R-r, 



4^ + 4/?: 



cos (261 - 26ij 



T4 



cos (20 - 36ij + 0k) ^cos (261 - 30^ + 0jk) 
\R-r,\'\R-rk\ ~ |il - r.^r, - r,.| 

cos(26ij - 26'fe) cos(46i-26'j -26*^,) 



\R-r,f\R-rkf 



\R-r,f\R^rkf 



(111) 



sin(26'- 36'j) 



r2 /l^ I r> |3 



y = y ± 



fcC 7 7 —I— C 



, 2r 



AT 

E 



cos(26i-36ij) 



Ii2-r,, 



2 sin (26* - 30J 

IV 



4 2r^ . cos (26* - 30i 

= y^T ^ — J — 73- 

iR-rA 



(112) 



VII. NUMERICAL RESULTS 



We propose to use the above-described transformation as the basis of a numerical method, whereby we make this 
transformation, integrate the reduced system, and invert the transformation. The overall efficacy of this Lie transform 
method will depend to a large degree on various factors such as the number, configuration, and circulation of vortices 
in the system. To probe the essential behavior of this method, we will consider the simplest case, that of three vortices 
of equal circulation, two of which are much closer to one another. This enables us to unambiguously state that our 
small perturbation parameter is e = |r|/|ii — rj|. For reference, we used a Runge-Kutta-Fehlberg (4, 5) adaptive 
time-stepping numerical integrator. 



A. Efficiency 



We would like to measure how much faster an algorithm based on the dynamics Eq. ( 106 ) and transformations 
Eq. (110 - 112) would run when compared to one based on the original point-vortex dynamics Eq. ([T]). The pertinent 
quantities we want to track are the number of calls each method makes to the base adaptive timer stepper. In 
particular, we will consider the ratio of these two numbers, i.e., the speedup factor, and how it behaves when e 
becomes small and our transformations therefore become more applicable. Since the time that it takes a vortex in the 
dimer to rotate through a given angle scales as the dimer separation squared Eq. (54 55|, we would roughly expect 



the speedup factor to scale as for all orders of Lie transforms. As one can see in Figure [ij this is very nearly the 
case. The magnitude of the speedup factor will certainly depend on the tolerance that one sets. However, increasing 
or decreasing the tolerance will not affect how the speedup factor scales with e. Neither will this scaling be modified 
by the contributions to the Lie transform algorithm run-time from the transformations themselves. While these do 
involve pairwise sums, they are calculated only when saving the state of the system is needed, and therefore don't 
constitute anything worse than a small multiplicative contribution to the overall run-time. Indeed their inclusion is 
only an additive contribution to the run-time if you are interested only in initial and final vortex configurations. 
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Speedup and Efficiency 




0.05 0.1 0.15 0.2 0.25 0.3 



FIG. 1. Tlie vertical axis shows the ratio of the number of calls to the adaptive time-stepper for the Regular Method (Nr) to 
the Lie Method {Nl), i.e. the speed-up factor. This is plotted against e, the vortex dimer separation |r„ — rm\ divided by the 
separation between the remaining vortex and the dimer \R — rj \ . A log- log version of this graph is inset. This corresponds 
to a power law -j^ oc t" where a « —1.8 



B. Accuracy 

We characterize the accuracy of the Lie transform algorithm by comparing how two identical point-vortex systems 
differ after evolving one with the regular method, and the other with the Lie method. Their deviation is quantified by 
the natural norm of the difference of their positions in configuration space. To have the regular method be a stand-in 
for the actual motion we must set the tolerance to be very small, which forces us to consider only the deviations 
that occur at small times. These deviations scale linearly with time, so the quantity of interest is really the deviation 
divided by time. Furthermore, there is a natural separation of the overall time-normalized deviation into that of 
the positions of the dimer vortices with respect to their center of circulation and that of the positions of the dimer 
itself and extra vortex. The time- normalized deviation of the first group is plotted against e in figure [2]for three 
orders of the Lie transformation, while that of the second group is plotted in a similar manner in figure ^ though 
for only two orders of the Lie transformation. Notice that because the constituents of the vortex pair have equal 
circulation, there is no 3rd order transformation to consider. The most important thing to note about each graph 
is that the time-normalized deviation scales as e", where a increases with the order of Lie transform used. That is, 
the Lie method gets more accurate as the dimer pair separation decreases and as we include higher orders of the Lie 
transform. 



C. Energy Conservation 



With respect to the close vortices circulating in a dimer, most numerical integration methods will systematically 
accrue errors that increase their separation. This is equivalent to decreasing the total energy, which is forbidden in an 
autonomous Hamiltonian system Eq. (16 1. With the regular integration method, the decrease in total energy per unit 
time is roughly proportional to the specified tolerance. More importantly, we can see from Figure |4] that this change 
in total energy per unit time roughly sales as e~^. This indicates that the non-constancy of the energy will certainly 
become a problem at small separations, no mater the tolerance. The performance of our Lie method is in marked 
contrast to that of the regular method. At every order of the Lie transform, the average energy remains constant. 
On top of this desired constant energy there are sinusoidal oscillations with a frequency commensurate with that of 
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Accuracy of the Dimer Internal Motion 




0.05 0.1 0.15 0.2 0.25 



e 



FIG. 2. The vertical axis shows the natural norm of the difference between the final configuration space points evolved under 
the Lie and Regular methods, normalized by the time elapsed. In particular, in this graph we consider the positions of only 
the dimer vortices with respect to their center of circulation. The three lines, from top to bottom, represent the Lie transform 
at 0th, 2nd, and 4th orders respectively. The inset log-log plot shows that they all follow power laws with exponents 1, 3, and 
5 respectively. 



Accuracv of the Rest of the System 



0.003 



0.0025 



0.002 



0.0015 



0.001 



0.0005 




1 y 




0.25 



FIG. 3. The vertical axis shows the natural norm of the difference between the final configuration space points evolved under 
the Lie and Regular methods, normalized by the time elapsed. In particular, in this graph we consider the positions of only the 

dimer and other vortex. The two lines, from top to bottom, represent the Lie transform at 0th - 2nd and 4th orders respectively. 
The inset log-log plot shows that they both follow power laws with exponents 2 and 4 respectively. 
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the dimer's rotation and an amplitude that depends on the order of the Lie transform and e. As indicated in Figure 
|4]the amphtude scales a,s A (x e". For the Lie transform at orders 0, 2, and 4, we have that a is nearly 2, 4, and 6 
respectively. Thus the accuracy of energy conservation under each order of the Lie transform method becomes better 
as e decreases. 



Energy Conservation 
0.01 F < < < < 




le-10 r 



le-11 I < < < < < 1 

0.02 0.04 0.06 0.08 0.1 0.2 0.3 

e 

FIG. 4. The vertical axis of this log-log plot shows two different items. The first is the change in energy per unit time, 
for evolution with the regular method. This corresponds to the line which increases as e decreases. Note that this is a power 
law, ^ oc e", with a of roughly -3. The second item is the amplitude. A, of oscillations in the energy about the fixed average 
for evolution with the Lie method. This corresponds to the three line which decrease as e decreases. Note that each follows a 
power law A oc e°'. From top to bottom, which corresponds to the Lie method at 0th, 2nd, and 4th orders, a is roughly 2, 4, 
and 6 respectively. 



D. Algorithmic Concerns 



In our implementation of this algorithm, the Lie transform method is triggered only when e dips below a pre- 
determined value. One then transforms the variables, evolves the new Hamiltonian, and waits for e to go above 
some de-triggering value before transforming back. It is fruitful to implement this algorithm in an object oriented 
manner such that the collection of N vortices is itself an object. When the method is triggered, we simply make a 
new (N-l)-vortex object with the transformed variables, and keep track of J and 9 externally. This view needs to 
be modified when you would like to use the transformation at fourth order or higher, where the modification to the 
Hamiltonian forces us to integrate the internal vortex dimer variables in conjunction with all the other variables. If 
using the Lie transforms only to third order does not seem to be much of a sacrifice, then the N-vortex object allows 
us to deal with more than one dimer pair at the same time. In particular we would just "pop" down one more level 
to an (N-2)-vortex system where the integrable dimer motions are treated separately. In particular, this would allow 
us to treat our simplest case of three vortices in a completely integrable fashion, where after considering the close 
vortices as a dimer, we could then look at the dimer and remaining vortex as a dimer itself. In this way, our method 
finds subsections of the phase-space that are integrable in a particular manner and exploits that for numeric gain. So, 
one must choose at third order whether one wants increased accuracy or the ability to treat other vortex pairs in a 
similar manner. 
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VIII. SUMMARY AND CONCLUSIONS 



The 2D poiiit-vortcx model is a very useful, albeit simple, approximation to many important physical systems. 
The forward time dynamics are easily obtained from a coupled nonlinear system of ODEs. However several problems 
arise when using an adaptive time-stepping method to numerically integrate the system while two like-signed vortices 
are very close compared to other vortex separations. If one keeps a small tolerance, the time-step will become pro- 
hibitively small and unnecessarily slow down the integration of the whole system. On the other hand, if one increases 
the tolerance to alleviate this problem, not only does the accuracy suffer, but the overall energy of the supposedly 
Hamiltonian system tends to drift. 

To address these issues, we transformed the system such that it separated into two separate components: the 
motion of the vortex pair, or dimer, about its center of circulation and the motion of the dimer itself, considered as a 
single vortex, along with the rest of the vortices. Wc viewed this as an (N-l)-vortcx system with the same dynamics, 
where one of the vortices, the dimer, has an internal "spin" degree of freedom. The corrections to this appealing 
picture were fortunately ordered in a small parameter. We then used Lie transform perturbation theory to develop a 
near identity canonical transformation that preserved this picture order by order in the small parameter. At fourth 
order, a slight alteration to the dynamics was necessary, though this affected only the dimer's spin degree of freedom 
by slightly lowering its rotation rate. With this inclusion, the (N-l)-vortex picture was altered only in so far as the 
spin of the dimer weakly coupled to a new field produced by the (N-1) vortices. 

When quantifying the efficacy of the algorithmic implementation of our transformations, the Lie method, we chose 

the particularly simple example of three vortices of equal circulation with one close pair. The first metric we con- 
sidered was the speed-up factor, which measures the ratio of the number of calls to the numerical integrator made 
by the regular method versus those made by the Lie method. Most importantly, we found this value to scale as 
e~i-^, which indicates that the regular method run-time diverges as e becomes small, since the Lie method run-time 
remains constant. Next we considered the accuracy of the Lie method versus the real motion as given by the regular 
method with a very low tolerance. The accuracy is quantified by the distance between two initially identical phase 
space points after evolving under the regular and Lie methods for a unit time. This was further broken down into 
the the phase space variables of the dimer relative to the center of circulation, and of the center of circulation and 
remaining vortex. For both cases wc found that the measure of accuracy gets smaller, i.e. more accurate, as a e 
decreases. As expected, the relationship is a power-law, with the integer power depending on the Lie transform 
order that is used. Additionally we examined the total energy of the system. The Lie method has zero drift in 
the value of the Hamiltonian as well as very small amplitude oscillations. These oscillations drastically decrease in 
amplitude as e decreases, and do so in a manner that gets better with inclusion of successively higher orders of the Lie 
transformation. This is to be contrasted with the consistent decrease in energy that occurs with the regular method 
when using an insufficiently low tolerance. Finally, we described how this transformation could be used as the basis 
of an object-oriented method for accurate numerical integration of point-vortex dynamics. 
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